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Abstract 

We report on a program for the numerical evaluation of divergent multi-loop integrals. The 
program is based on iterated sector decomposition. We improve the original algorithm of 
Binoth and Heinrich such that the program is guaranteed to terminate. The program can be 
used to compute numerically the Laurent expansion of divergent multi-loop integrals regu- 
lated by dimensional regularisation. The symbolic and the numerical steps of the algorithm 
are combined into one program. 



PROGRAM SUMMARY 



Title of program: sector_decomposition 
Version: 1.0 
Catalogue number: 

Program obtained from: http : / /wwwthep .physik .uni-mainz . de/~stefanw/ software .html 
E-mail: stef anw@thep .physik . uni-mainz . de, bogner@t hep .physik .uni-mainz . de 
License: GNU Public License 
Computers: all 
Operating system: Unix 
Program language: C++ 

Memory required to execute: Depending on the complexity of the problem. 

Other programs called: GiNaC, available from http : / / www . ginac . de, 
GNU scientific library, available from http : //www. gnu . org/software/gsl. 

External files needed: none 

Keywords: Multi-loop integrals. 

Nature of the physical problem: Computation of divergent multi-loop integrals. 
Method of solution: Sector decomposition. 

Restrictions on complexity of the problem: Only limited by the available memory and CPU time. 
Typical running time: Depending on the complexity of the problem. 
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1 Introduction 



The calculation of higher-order corrections in perturbation theory in particle physics relies to a 
large extent on our abilities to compute Feynman loop integrals. Any calculation of loop integrals 
is complicated by the fact, that these integrals may contain divergences, ultraviolet or infrared in 
nature. Dimensional regularisation is usually employed to regulate the divergences and one seeks 
the Laurent expansion of Feynman integrals in D = 4 — 2e dimensions in the small parameter £. 

Automated numerical algorithms, which permit the computation of the coefficients of the 
Laurent series for certain phase space points are of great use as an independent check of ana- 
lytical calculations. There are several methods available for this task. All methods have to face 
the challenge to disentangle overlapping singularities. In this paper we focus on sector decom- 
position [1-4]. Other possibilities include a method based on difference equations [5-8] or the 
evaluation of Mellin-Barnes integrals [9, 10]. 

Sector decomposition is a convenient tool to disentangle overlapping singularities. Binoth 
and Heinrich [3,4] gave an algorithm for the automated computation of the Laurent series of 
multi-loop integrals. The method has recently also been applied to phase space integrals [11-13] 
and complete loop amplitudes [14, 15]. 

However, the original formulation of Binoth and Heinrich has a drawback: The disentan- 
glement of overlapping singularities is achieved through recursive sector decomposition. If this 
process terminates, the algorithm gives the correct answer. But there is no guarantee that the 
algorithm stops after a finite number of steps. In fact, one can find relatively simple counter- 
examples, where the original algorithm leads to an infinite loop. In this paper we improve the 
method and present an algorithm, which is guaranteed to terminate. This is possible, since we 
can relate our problem to a well-known problem in algebraic geometry: The resolution of sin- 
gularities of an algebraic variety over a field of characteristic zero by a sequence of blow-ups. 
In view of [16] a connection between Feynman integrals and algebraic geometry comes to no 
surprise. The mathematical problem can be solved, which was first shown by Hironaka [17]. Al- 
though the proof by Hironaka is non-constructive and therefore cannot be used as an algorithm, 
in subsequent years constructive methods have been developed [18-28]. Let us mention for com- 
pleteness that the mathematical problem is more general, the disentanglement of singularities of 
multi-loop integrals corresponds to the sub-case, where the algebraic variety is generated through 
a single polynomial. We have implemented three different strategies for the resolution of the sin- 
gularities. All of them are guaranteed to terminate and yield the correct answer. They differ in 
the required CPU-time. 

The algorithm for sector decomposition can be divided into two stages. The first stage in- 
volves symbolic computer algebra, while the second stage is based on a numerical integration. 
Partially due to this "double nature" of the algorithm, no public program is available up to now. 
Given the importance and the popularity of the algorithm, this gap should be closed. With this 
article we provide such a program. 

A user-friendly implementation of the algorithm would require that the different aspects of 
the method, e.g. the symbolic and the numerical computations, are treated within the same frame- 
work. The library GiNaC [29] allows the symbolic manipulation of expressions within the pro- 
gramming language C++. It offers the tools to combine symbolic and numerical computations 
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in one program. Furthermore there exist specialised libraries for an efficient numerical Monte 
Carlo integration [30-32]. Our implementation is therefore in C++, using the libraries GiNaC 
for the symbolic part and the GNU Scientific Library for the numerical part. All stages of the 
algorithm are implemented within one program. 

This article is organised as follows: In section 2 we review basic facts about multi-loop 
Feynman integrals. Section 3 gives a brief description of the algorithm for sector decomposition. 
In section 4 we discuss in detail strategies, which guarantee a termination of the iterated sector 
decomposition. Section 5 gives an overview of the design of the program while section 6 is of a 
more practical nature and describes how to install and use the program. This section provides an 
example and compares the efficiency of the different strategies. Finally, a summary is provided in 
section 7. In an appendix we describe technical details of the implementation and provide a proof 
for one of the strategies (the proofs for the two other strategies can be found in the literature). 



2 Feynman integrals 

In this section we briefly recall some basic facts about multi-loop Feynman integrals. A generic 
scalar /-loop integral Ig in D dimensions with n propagators corresponding to a graph G is given 

by 

r 1 d D k n 1 

ig = n-4n ( 2? 2,, (d 

J r=i in? j=\ (-qj + mjyj 

The momenta qj of the propagators are linear combinations of the external momenta and the loop 
momenta. After Feynman parameterisation one arrives at 

nr(v 7 ) io <=i V=i J 7 

7=1 

with v = Vi + ... + V n . The functions 11 and J depend on the Feynman parameters. If one 
expresses 

j^Xji-qj + mj) = - j^j^krMrsks+j^lkr-Qr-J, (3) 

7=1 r=ls=\ r=\ 

where M is a / x I matrix with scalar entries and Q is a /-vector with four-vectors as entries, one 
obtains 

U = det(M), 7 = det(M) (-J + QM- l Q) . (4) 

Alternatively, the functions 11 and J can be derived from the topology of the corresponding 
Feynman graph G. Cutting / lines of a given connected /-loop graph such that it becomes a 
connected tree graph T defines a chord c{T,G) as being the set of lines not belonging to this 
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tree. The Feynman parameters associated with each chord define a monomial of degree /. The 
set of all such trees (or 1 -trees) is denoted by 1\. The 1 -trees T e 1\ define U as being the sum 
over all monomials corresponding to the chords c(T,G). Cutting one more line of a 1-tree leads 
to two disconnected trees (T\,T2), or a 2-tree. T2 is the set of all such pairs. The corresponding 
chords define monomials of degree I + 1 . Each 2-tree of a graph corresponds to a cut defined 
by cutting the lines which connected the two now disconnected trees in the original graph. The 
square of the sum of momenta through the cut lines of one of the two disconnected trees T\ or T2 
defines a Lorentz invariant 



s T = 



Pj 



(5) 



Kjec(T,G) 



The function 7o is the sum over all such monomials times minus the corresponding invariant. 
The function f is then given by 7o plus an additional piece involving the internal masses rrij. In 
summary, the functions 11 and J are obtained from the graph as follows: 



U = 



H = 



7 = 



e n 

Te?x j€C{T,G) 



{T u T 2 )eT 2 jec{T u G) 
n 

7o + uY, Xjm) . 
7=1 



EI X J (-^1) , 



(6) 



In general, 11 is a positive semi-definite function. Its vanishing is related to the UV sub- 
divergences of the graph. Overall UV divergences, if present, will always be contained in the 
prefactor T(v — ID/ 2). In the Euclidean region, f is also a positive semi-definite function of the 
Feynman parameters xj. The vanishing of f is related to infrared divergences. Note that this is 
only a necessary but not sufficient condition for the occurrence of an infrared singularity. If an 
infrared singularity occurs or not depends in addition on the external kinematics. 



3 A review of the algorithm for sector decomposition 

In this section we briefly describe the original algorithm of Binoth and Heinrich for iterated 
sector decomposition. A detailed discussion on strategies for choosing the sub-sectors is deferred 
to section 4. The program calculates the Laurent expansion in £ of integrals of the form 

/ d-*8(i-2>) (f\xr £b ) m p M Cj+£di - ( 7 ) 

x/>o i=1 v=i /;=i 

The integration is over the standard simplex. The a's, b's, c's and d's are integers. The P's are 
polynomials in the variables x\, x n . The polynomials are required to be non-zero inside the 
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integration region, but may vanish on the boundaries of the integration region. Note that the 
program allows a product of several polynomials and that it is not required that the polynomials 
are homogeneous. 



Step 0: We first convert all polynomials to homogeneous polynomials. Due to the presence 
of the delta-function we have 

1 = X1+X2 + ...+x n . (8) 

For each polynomial Pj we determine the highest degree hj and multiply all terms with a lower 
degree by an appropriate power of x\ +X2 + ... +x n . As a result, the polynomial Pj is then ho- 
mogeneous of degree hj. 



Step 1 : We decompose the integral into n primary sectors as in 



f d n x = £ f d n x n e(jc/>jc,-). 
{„ i=i in i=\,m 



(9) 



xj>0 xj>0 



#1 

In the Z-th primary sector we make the substitution 

Xj = xix'j for j 7^ /. (10) 
As each polynomial Pj is now homogeneous of degree hj we arrive at 

/ d»x8(i-j>)( fl e(xi>*)) (fl-? +eh )fl[PM Cj+£dj = 

in i=i \i=h¥l J V=i ) i=\ 



Xj >0 



[( fl dx iX ? +e A(l+ £ x^\ f\[Pj(x l ,...,Xj- h l,x j+h ...,x n )} Cj+edj ,(n) 

o v=h¥i J V j=hm J 7=1 



where 



c = -n - £ (o ; + efc,-) - £ /z j (c 7 - + edj) . (1 2) 

(=1 7=1 

Each primary sector is now a (n — 1) -dimensional integral over the unit hyper-cube. Note that in 
the general case this decomposition introduces an additional polynomial factor 

1+ £ xi] . (13) 

However for Feynman integrals of the form as in eq. (2) we always have c = and this factor 
is absent. In any case, this factor is just an additional polynomial. In general, we therefore deal 
with integrals of the form 

l 

f d n xf[x? +£b ' fl [Pj(x)} Cj+£dj . (14) 
o !=1 j =1 



Step 2: We decompose the primary sectors iteratively into sub-sectors until each of the polyno- 
mials is of the form 

P = x™K..x%»(l+P f (x)), (15) 

where P'(x) is a polynomial in the variables xj. If P is of the form (15), we say that P is 
monomialised. In this case the monomial prefactor x™ 1 ...x% n can be factored out and the re- 
mainder contains a constant term. To convert P into the form (15) we choose a subset S = 
{oci, OCfc} C {1, ...n} according to one of the strategies outlined in section 4. We decompose 
the ^-dimensional hypercube into k sub-sectors according to 




In the Z-th sub-sector we make for each element of S the substitution 

*ctj = *a/ a , forzV/. (17) 

This yields a Jacobian factor x^ 1 . This procedure is iterated, until all polynomials are of the 
form (15). The strategies discussed in section 4 guarantee that the recursion terminates. At the 
end all polynomials contain a constant term. Each sub-sector integral is of the form as in eq. (14), 
where every Pj is now different from zero in the whole integration domain. Hence the singular 
behaviour of the integral depends on the a,- and bi, the a, being integers. 



Step 3: For every xj with aj < we have to perform a Taylor expansion around xj = in order to 
extract the possible £-poles. If we consider for the moment only one parameter xj we can write 
the corresponding integral as 







where we defined 1^ = d/dxjl(xj) 



. The remainder term 



(18) 



xj=0 



i( R \ Xj ) = i(xj)- £ (19) 

P =o P- 

does not lead to £-poles in the Xj-integration. The integration in the pole part can be carried out 
analytically: 

fdxjx'^illP) = -i (20) 

J J J pi dj + bjE + p + l p\ 
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This procedure is repeated for all variables Xj for which aj < 0. 

Step 4: All remaining integrals are now by construction finite. We can now expand all ex- 
pressions in a Laurent series in £ 

B 

£c i e i + 0(e B ) (21) 

i=A 

and truncate to the desired order. 

Step 5: It remains to compute the coefficients Q. These coefficients contain finite integrals, 
which we evaluate numerically by Monte Carlo integration. 

4 Strategies for choosing the sub-sectors 

In this section we discuss strategies for choosing the subset S which defines the sub-sectors, such 
that the iterated sector decomposition terminates. Binoth and Heinrich proposed to determine 
a minimal subset S = {oci, oc^} such that at least one polynomial Pj vanishes for x ai = ... = 
*a k = 0. By a minimal set we mean a set which does not contain a proper subset having this 
property. If S contains only one element, S = {a}, then the corresponding Feynman parameter 
factorises from Pj. A relative simple example shows, that this procedure may lead to an infinite 
loop: If one considers the polynomial 

P(xi,X2,X3) =X\xl+X2+X2X3, (22) 

then the subset S = {1, 2} is an allowed choice, as P(x\ = 0, X2 = 0, x$) = and S is minimal. 
In the first sector the substitution (17) reads x\ = x\ = x' 3 . It yields 

D / \ / i2 i 12 12 , / / / I ( 12 , I 12 , I i \ / r> I i i i\ 

r \X\ , X2j Xt, ) — -^1-^3 \X^X2 ~rJCjJC2-^3 — -^1 1-^3 _ ' _ '^2"^3/ — X^r IXj , X3, X2J • 

The original polynomial has been reproduced, which leads to an infinite recursion. The problem 
of a non-terminating iteration has already been noted in ref. [4]. To avoid this situation we need 
a strategy for choosing S, for which we can show that the recursion always terminates. This is a 
highly non-trivial problem. It is closely related to the resolution of singularities of an algebraic 
variety over a field of characteristic zero. Fortunately, mathematicians have found a solution for 
the latter problem and we can adapt the mathematical solutions to our problem. We will not go 
into any details of algebraic geometry. Instead we present in section 4.1 the polyhedra game, 
introduced by Hironaka to illustrate the problem of resolution of singularities. The polyhedra 
game can be stated with little mathematics. Any solution to the polyhedra game will correspond 
to a strategy for choosing the subsets S. In sections 4.2 to 4.4 we present three winning strate- 
gies for the polyhedra game. All three strategies ensure that the iterated sector decomposition 
terminates and lead to the correct result. However, the number of generated sub-sectors (and 
therefore the efficiency of the method) will vary among the different strategies. Common to all 
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strategies is a sequence of positive numbers associated to the polynomials. All strategies enforce 
this sequence to decrease with each step in the iteration with respect to lexicographical ordering. 
As the sequence cannot decrease forever, the algorithm is guaranteed to terminate. The actual 
construction of this sequence will differ for different strategies. 

The lexicographical ordering is defined as follows: A sequence of real numbers (ai, ...,a r ) is 
less than a sequence (b\, ...,b s ) if and only if there exists a k e N such that aj = bj for all j < k 
and a k < b^, where we use the convention that aj = for j > r. 



4.1 Hironaka's polyhedra game 

Hironaka's polyhedra game is played by two players, A and B. They are given a finite set M 
of points m = (mi, m n ) in N" , the first quadrant of W. We denote by A C W]_ the positive 
convex hull of the set M. It is given by the convex hull of the set 

U (m + R n + ). (24) 

The two players compete in the following game: 

1. Player A chooses a non-empty subset S C {1, n}. 

2. Player B chooses one element i out of this subset S. 

Then, according to the choices of the players, the components of all (mi , m n ) G M are replaced 
by new points (m[, ...,m' n ), given by: 

jf = m h if iV z \ 

i = E m J- c ' ( 25 ) 



m 
m 



where for the moment we set c = 1. This defines the set M' . One then sets M = M' and goes 
back to step 1. Player A wins the game if, after a finite number of moves, the polyhedron A is of 
the form 

A = m + R n + , (26) 

i.e. generated by one point. If this never occurs, player B has won. The challenge of the polyhe- 
dra game is to show that player A always has a winning strategy, no matter how player B chooses 
his moves. 

Let us discuss the relation of Hironaka's polyhedra game to the sector decomposition of 
multi-loop integrals. Without loss of generality we can assume that we have just one polynomial 
P in eq. (7). (If there are several polynomials, we obtain a single polynomial by multiplying them 
together. As only the zero-sets of the polynomials are relevant, the exponents can be neglected.) 
The polynomial P has the form 

p OT (0 OT (0 (<) 
(=1 
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The n-tuple m (0 = (m { -\ ...,mi°J defines a point in N" andM= jmW, ...m^j is the set of all 
such points. Substituting the parameters xj according to equation (17) and factoring out a term 
yields the same polynomial as replacing the powers mj according to equation (25). In this sense, 
one iteration of the sector decomposition corresponds to one move in Hironaka's game. Reducing 
P to the form (15) is equivalent to achieving (26) in the polyhedra game. Finding a strategy which 
guarantees termination of the iterated sector decomposition corresponds to a winning strategy for 
player A in the polyhedra game. Note that we really need a strategy that guarantees player A's 
victory for every choice player B can take, because the sector decomposition has to be carried 
out in every appearing sector. In other words, we sample over all possible decisions of B. 



4.2 Strategy A 

Zeillinger's strategy [27, 28] for player A is conceptionally the simplest. A drawback of this 
strategy is that it leads to a large number of sub-sectors and is therefore quite inefficient. Denote 
by Mq the set of corners of A, Mq C M. Define W to be the set of vectors connecting these 
corners: 

W = | v = m W - ra W : ra W , ra W e M c } . (28) 

We define the length L(v) of the vector v <EW by 

L(v) = max v; — min v,-. (29) 

\<i<n \<i<n 

We define the multiplicity N(v) of the vector v G W by the number of its components equal to its 
minimal or maximal component: 



N(v) = # < j : 1 < j < n,vj = min v,- or vj = max v,- > . 

1 l<i<n \<i<n J 



(30) 



From the set W we single out one vector w for which the sequence of the two numbers L (w) and 
Af (w) is minimal with respect to lexicographical ordering: 

(L(w),N(w)) < lex (L(v),7V(v)), VveW. (31) 

Player A chooses then S = {k, I}, where k and / are defined by 

Wk = min Wi and w/ = max W{. (32) 

\<i<n \<i<n 

This choice reduces the sequence of the three numbers 

^ _ f (0, 0, 0) if A is monomial, 

~ \ (#M c ,L(v M ), AT(vm)) otherwise, ! J ' ; 

with respect to lexicographical ordering. #Mc denotes the number of corners of A. A proof that 
this strategy is a solution to Hironaka's polyhedra game can be found in [27, 28]. 
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4.3 Strategy B 

Spivakovsky's strategy [18] was the first solution to Hironaka's polyhedra game. To state the 
algorithm, we need a few auxiliary definitions: We define co (A) e R + as the vector given by the 
minima of the individual coordinates of elements in A: 

co, = min{v ; 1 v e A} , i = 1, n. (34) 

Furthermore we write A = A — co (A) and V/ = V; — CO/. For a subset r C {1, ...,«} we define 

d r {A) = minj|>;|veAj and d{A) = d {l ^ n} (A) . (35) 

We then define a sequence of sets 

(7o,A ,/i,Ai,...,7 r ,A r ) (36) 

starting from 

7 = {1, ...,«}, A =A. (37) 

For each A^ we define a set Hk by 

74 = < jel k \ 3 v G A fc such that £v/ = d(A fc ) andv^ol. (38) 
I ie/* J 

4 + i is given by 

4+i = 4\74- (39) 

In order to define A^+i we first define for the two complementary subsets 74 and 7^+1 of 7^ the 
set 



M Hk = ivGR* | £ v / < 1 1 



(40) 



and the projection 



T 3 ^ : 



^(«»P) = TZ1ri' PGR?, IPI=EPy. (41) 

Then A fc+1 is given by 

4t+i = P^^n^^-yUA^j, (42) 

where A^ = A^ — co(A^). The sequence in eq. (36) stops if either d (A r ) = or A r = 0. Based on 
the sequence in eq. (36) player A chooses now the set 5 as follows: 
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1. If A r = 0, player A chooses S = {1, n}\I r . 

2. If A r 7^ 0, player A first chooses a minimal subset T r C 7 r , such that Ejer r v j > 1 f° r ai l 
V G A r and sets 5 = ({1, ...,n}\/ r )ur r . 

To each stage of the game (i. e. to each A), we can associate a sequence of 2r + 2 numbers 

8(A) = (d(A),#h,d(A 1 ),...,#I r ,d(A r ),d(A r )), (43) 

adopting the conventions = and d (0) = °°. The above strategy forces 8 (A) to decrease with 
each move in the sense of a lexicographical ordering. Further, it can be shown that 8 (A) cannot 
decrease forever. Hence player A is guaranteed to win. The proof is given in [18]. 



4.4 Strategy C 

In this section we present a strategy inspired by the proof of Encinas and Hauser [24, 26]. We 
keep as much as possible the notation of the previous section. Similar to the previous section we 
define a sequence 

(c_i,7o,A ,co,/i,Ai,...,c r _i,7 r ,A r ) (44) 

starting from 

c_i=0, 7 = {1, ..,«}, A = A (45) 

Compared to eq. (36) we added the numbers c_i,co, ...,c r _i to the sequence. Again we define 
H k by 



(46) 



H k = \jel k \ 3veA k such that £ v,- = d (A k ) and v, ^ > . 
We set Ck = d (A^) and define for a set A^ a companion set A^ by 

y A k otherwise. 

Let T k C H k be a subset of H k . The subset T k is chosen according to a rule 7?: 7^ = R(H k ). The 
only requirement on the rule R is, that is has to be deterministic: If 

T = R(H) and T 1 = R(H') and 77 = 77 ; =>• 7 1 = 7'. (48) 

We define 

A Tk = c k P Tk [M Tk n^-). (49) 



Ck 
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We then set 



4+i = h\Tk, 

A k+ i = A Tk . (50) 

The sequence in eq. (44) stops if either d (A r ) = or A r = 0. Based on the sequence in eq. (36) 
player A chooses now the set S as follows: 

1. If A r = 0, player A chooses S = {1, n}\I r . 

2. If A r 7^ 0, player A first chooses a minimal subset T r C I r , such that Ejer r v j > c f--i for all 
VGA r and sets S = ({1, «} \/ r ) UT r . 

This strategy reduces the sequence 

/(A) = (d(Ao),{M -#H ),d(A 1 ),...,{#I r - l -#H r -l),d(A r ),d{A r )) (51) 

with respect to lexicographical ordering. The strategy C is similar to strategy B, the major dif- 
ferences are in the choice of the companion set eq. (47) as compared to eq. (42) and the freedom 
to choose a subset 7^ C H K instead of H^. We provide a proof for this strategy in appendix B. 

5 A description of the program 

The program is written in C++. The main routine to compute an integral of the form as in eq. (7) 
is the function do_sector_decomposition. The arguments are as follows: 

monte_carlo_result do_sector_decomposition ( 

const integration_data & global_data, 

const integrand & integrand_in, 

const monte_carlo_parameters & mc_parameters, 

int verbose_level = 0); 

The input are three structures, integration_data, integrand, monte_carlo_parameters, 
which are described in detail below, and an optional parameter verbose_level. With the help of 
the optional parameter verbose_level one can choose the amount of information the program 
prints out during the run. The function returns a structure monte_carlo_result, which again 
is described below. The strategy for the iterated sector decomposition is selected by the global 
variable CHOICE_STRATEGY. The keywords for the available strategies are 

STRATEGY_A, STRATEGY_B, STRATEGY_C, STRATEGY_X . 

The first three strategies are described in detail in the previous section and are guaranteed to 
terminate. The last one is a heuristic strategy, for which we neither have a proof that it terminates, 
nor do we know a counter-example which would lead to an infinite recursion. This strategy 
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chooses the smallest set S such that the maximal power of a Feynman parameter can be factored 
out in each sub-sector. It is included for the following reason: If this strategy terminates, the 
number of the generated sub-sectors tends to be smaller than the corresponding numbers for the 
other strategies. The default is 

CHOICE_STRATEGY = STRATEGY_C; 

The class integration_data contains the data which will not be modified by the algorithm for 
the sector decomposition. It has a constructor of the form 

integration_data (const std : : vector<GiNaC : : ex> & list_feynman_parameter, 
GiNaC::ex epsilon, int order); 

where list_f eynman_parameter is a vector holding the symbols of the Feynman parameters, 
epsilon is a symbol corresponding to e in eq. (7) and order defines which term of the Laurent 
series in £ should be computed. 

The integrand 



V=i / j =i 

is encoded in the class integrand. This class has a constructor of the form 

integrand (const std : : vector<exponent> & nu, 

const std: :vector<GiNaC: :ex> & poly_list, 
const std: :vector<exponent> & c) ; 

where nu is a vector of size n holding the exponents <2 ; - + £& ; -. poly_list is a vector of size r, 
holding the polynomials Pj in the Feynman parameters. The corresponding exponents are stored 
in the vector c, again of size r. As exponents are generally of the form a + be, a special class is 
available for them: 

exponent (int a, int b) ; 

In applications one encounters often integrands where the explicit powers of the Feynman pa- 
rameters are missing. For integrands of the form 




(52) 



r 



Cj+zdj 



n [ p M 



(53) 



7=1 



there is a simpler constructor of the form 



integrand(size_t n, const std : : vector<GiNaC : : ex> & poly_list, 
const std: :vector<exponent> & c) ; 
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Here n is the number of Feynman parameters. 

Parameters associated to the Monte Carlo integration are specified with the help of the class 
monte_carlo_parameters. This class is constructed as follows: 

monte_carlo_parameters (size_t iterations_low, size_t iterations_high, 

size_t calls_low, size_t calls_high) ; 

The program uses the Vegas-algorithm [33, 34] based on an adaptive grid. The program does 
first a Monte Carlo integration with iterations_low iterations with calls_low function eval- 
uations each. This phase is solely used to adapt the grid. The numerical result of the Monte 
Carlo integration is then obtained from the second stage with iterations_high iterations of 
calls_high function calls each. 

The main function do_sector_de compos it ion returns the numerical results of the Monte Carlo 
integration in the class monte_carlo_result. This class has the methods 

class monte_carlo_result { 
public : 

double get_mean (void) const; 
double get_error (void) const; 
double get_chi_squared (void) const; 

}; 

which return the mean value of the Monte Carlo integration, the error estimate and the associated 
X 2 - 

6 How to use the program 

In this section we give indications how to install and use the program library. Compilation of 
the package will build a (shared) library. The user can then write his own programs, using the 
functions provided by the library by linking his executables against the library. 

6.1 Installation 

The program can be obtained from 

http : / /wwwthep .physik . uni-mainz . de/~stefanw/ software .html 

After unpacking, the library for sector decomposition is build by issuing the commands 

. / configure 
make 

make install 
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There are various options which can be passed to the configure script, an overview can be ob- 
tained with . /configure — help. 



After installation, the shell script sector_decomposition-conf ig can be used to determine 
the compiler and linker command line options required to compile and link a program with 
the library. For example, sector_decomposition-conf ig — cppflags will give the path to 
the header files of the library, whereas sector_decomposition-conf ig — libs prints out the 
flags necessary to link a program against the library. 

6.2 Writing programs using the library 

The following test program computes the first terms of the Laurent series of the massless double- 
box graph. The graph corresponds to the integral 



f = [X2X^(X4+X5 +X(, +Xj) + XsX(,(X[ + X% +X3 + X4) +X2X4X6 +X3X4X5] (s) + X\X4Xq( — t). 

The integral is computed for the point (s,t) = ( — 1,-1). 

♦include <iostream> 
♦include <vector> 
♦include <ginac/ginac . h> 

♦include " sect or_decomposit ion/ sector_decomposition .h" 

int main() 
{ 

using namespace GiNaC; 

using namespace sector_decomposition; 

CHOICE_STRATEGY = STRATEGY_X; 

symbol eps ( "eps" ) ; 
symbol s ("s") , t ("t") ; 

symbol xl ("xl") , x2("x2"), x3("x3"), x4("x4"), x5("x5"), x6("x6"), x7("x7"); 
const ex x[] = { xl, x2 , x3, x4 , x5, x6, x7 }; 

std: :vector<ex> parameters (7) ; 

for (int i=0; i<7; i++) parameters [i] =x [i] ; 




(54) 



with 



U = (xi+X2+X3)(x5+X(,+X-i) +X4(xi+X2+X3+X5+X(,+X-i), 



(55) 
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ex U = (xl+x2+x3) * (x5+x6+x7) + x4 * (xl+x2+x3+x5+x6+x7 ) ; 

ex F = (x2*x3* (x4+x5+x6+x7 ) +x5*x6* (xl+x2+x3+x4 ) +x2*x4*x6+x3*x4*x5) * (-s) 
+xl*x4*x7* (-t) ; 

F = F . subs (1st ( s == -1 , t == -1 )); 

std: :vector<ex> poly_list (2 ) ; 
poly_list [0] = U; 
poly_list [1] = F; 

std: :vector<exponent> c(2); 
c[0] = exponent ( 1, 3 ); 
c[l] = exponent ( -3, -2 ); 

integrand my_integrand (7 , poly_list, c) ; 

monte_carlo_parameters mc_parameters ( 5, 15, 10000, 100000 ); 

for (int order=-4; order<=0; order++) 
{ 

integration_data global_data (parameters, eps, order); 

monte_carlo_result res = 

do_sector_decomposition (global_data, my_integrand, mc_parameters) ; 

std::cout << "Order " << pow (eps, order) << ": " << res . get_mean ( ) 
<< " +/- " << res . get_error ( ) << std: :endl; 

} 

return 0; 

} 

Running the program will print out the following result: 

Order eps A (-4): 2.00001 +/- 9.25208e-05 
Order eps A (-3): -5.99992 +/- 0.000359897 
Order eps A (-2): -4.91623 +/- 0.00157598 
Order eps A (-l): 11.4958 +/- 0.00681643 
Order 1: 13.8236 +/- 0.0207286 

The running time is about 40 minutes on a standard PC. 
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Integral 


Strategy A 


Strategy B 


Strategy C 


Strategy X 


Bubble 


2 


2 


2 


2 


Triangle 


3 


3 


3 


3 


Tbubble 


58 


48 


48 


48 


Planar double-box 


755 


586 


586 


293 


Non-planar double-box 


1138 


698 


698 


395 



Table 1 : The number of generated sub- sectors for different loop integrals and different strategies. 

6.3 Further examples and performance 

The program comes with several examples: The one-loop two-point and three-point functions, as 
well as the following two-loop functions: The two-point function, the planar double-box and the 
non-planar double box. The corresponding Feynman diagrams for these examples are shown in 
fig. 1. For these examples the corresponding analytic results are known and can be found for the 




TBubble Planar double-box Non-planar double-box 



Figure 1: Feynman integrals included as examples in the distribution. 

two-loop integrals in [35-37]. We have used these results to verify the correctness of our code. 

The running time of the program is dominated by the numerical Monte-Carlo integration. 
The running time for the Monte-Carlo integration depends on the complexity of the integrand, 
which in turn is related to the number of generated sub-sectors. The number of generated sub- 
sectors is therefore a measure for the efficiency of the algorithm. Note that although different 
strategies lead to the same result for the Laurent expansion of multi-loop integrals, the number of 
generated sub-sectors may differ among the various strategies. In table 1 we compare the number 
of generated sub-sectors for the different strategies and for different loop integrals. We observe 
that strategy A performs worse than strategies B or C. Strategies A, B and C are guaranteed 
to terminate. Strategy X, for which we have no proof that the singularities are resolved after a 
finite number of iterations, terminates for the examples above and gives the lowest numbers for 
the generated sub-sectors. For this reason it is included in the program. A pragmatic approach 
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would try strategy X first. If the iterated sector decomposition terminates, one obtains the result 
efficiently. If not, one uses as fall-back option one of the strategies A, B or C. 

Going to a higher number of loops and more propagators, the memory requirements become 
an issue. However, the individual sub- sectors are independent and can be calculated once at a 
time. This reduces significantly the memory requirements. With this method we evaluated the 
massless on-shell triple box in about two days CPU time on a standard PC. 

6.4 Documentation 

The complete documentation of the program is inserted as comment lines in the source code. 
The documentation can be extracted from the sources with the help of the documentation sys- 
tem "doxygen" [38]. The program "doxygen" is freely available. Issuing in the top-level build 
directory for the library the commands 

doxygen Doxyfile 

will create a directory "reference" with the documentation in html and latex format. 

7 Summary 

In this article we considered the numerical computation of the coefficients of the Laurent ex- 
pansion of a divergent multi-loop integral, regulated by dimensional regularisation. The method 
is based on iterative blow-ups, also known as iterative sector decomposition. The algorithms 
we employed ensure that this recursive procedure terminates. We implemented the algorithms 
into a computer program, such that the symbolic and the numerical steps of the algorithms are 
contained in one program. 
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A Details on the implementation 

The numerical Monte Carlo integration is by far the most time-consuming part of the program 
and efficiency of the program in this part is of particular importance. 

The GiNaC-library provides a method for the numerical evaluation of a function, based on 
arbitrary precision arithmetic. For Monte Carlo integration, where a function needs to be evalu- 
ated many times, this is quite slow and therefore inefficient. It is also not needed, since statistical 
errors and not rounding errors tend to dominate the error of the final result. Therefore a different 
approach has been implemented for the numerical Monte Carlo integration: The function to be 
integrated is first written as C code to a file, this file is then compiled with a standard C compiler 
and the resulting executable is loaded dynamically (e.g. as a "plug-in") into the memory space 
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of the program and the Monte Carlo integration routine uses this compiled C function for the 
evaluations. 

The performance is further improved by optimising the C code of the function. This is done 
by replacing subexpressions, which occur syntactically more than once by temporary variables. 
In addition, subexpressions with a large number of operands are split into smaller pieces. This 
ensures that long expressions can be processed by the compiler. These optimising techniques are 
part of many commercial computer algebra systems, but are not part of the GiNaC-library. They 
have been added to this program. 



B Proof of strategy C 

In this appendix we provide a proof of strategy C. The structure of the proof is modelled on 
Spivakovsky's proof of strategy B. We keep the notations of section 4. 



Lemma 0: 

d(A k ) > c*_i. (56) 

Proof by induction: The claim is obvious for k = 0, since c_i =0. Assume that the hypothesis 
is true for k. With d (A k ) = d (A k ) + |co(A^)| we deduce that d (Ajj;) > c^. Using this it follows 
d (Aj k ) > Ck for all T k . Therefore d (A k+ i) > c k . 

Corollary: d (A^) = d (Ajj;) . For Ck > c k -\ there is nothing to prove. For c k < c^-i we have 
to consider the point c^(o(A fc ) /{ck-\ — c*). From |to(A^)| = d(A k ) — d (A^) > c k -\ — the 
claim follows. 



We call a set Sk permissible for Ak, if 

£v 7 - > c k -i forallveA^. (57) 

Lemma 1 : If S k +i is permissible for A^+i, then Sk = Sk+i U T k is permissible for A^ and ds k (A k ) = 
d (Ayt) . To prove this lemma we first show 

> c k -i-c k , (58) 

where (Oj are the components of (o(Ayt). For Ck>Ck-\ this is obvious, as the l.h.s. is non-negative. 

For c k < Ck-\ we have 

Ck (»(A k )eA c k . (59) 



Ck-i -c k 
Then 

Sk+i permissible for A^+i =>- ^ v ; — c k f° r all VGAj =^ ^ <°j — c k-\ — Ck- (60) 

jes k+1 un jes k 
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Therefore 



Sk+i permissible for A^+i 



jeS k+l i)T k 

£ v j > c k for all VEA k 

j£S k jeS k 

£ Vj > cn forallveAfc 

jes k 



(61) 



From the second step above 



£ V; > c k = d (A k ) for all V e A k 

jts k 



(62) 



it follows immediately that ds k (Ak) > ^ (A/t) . 



Corollary: The set S as defined in section 4.4 is permissible for A. For A r ^ the set T r is per- 
missible by construction. Then the repeated application of the above lemma shows that S is per- 
missible for A. For A r = it follows that ds r _j (A^_i) > c r _i and therefore ds r _j (A r _i) > c r _i. 
Therefore S r -i is permissible for A r _i and with the same argumentation as above it follows that 
S is permissible for A. 

Let now A C W]_. We define the transformation A' = o c s • (A) by 



v'j = Vj for j ^ i, 

Lemma 2: 

(a) d(A! k )<d(A k ). 

(b) If J (Ay = d (A*) then i £ ff* and C 

(c) If J (AQ = d (A k ) and fljt = H' k then 7* = T k \ 

(d) If J (A^) = d (Ayt) and 7^ = T k ' then the following diagram commutes: 



(63) 



h+1 



T e /t-i 



S *+l>' 



A' 



- A' 



- A' 



jfc+i 
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Proof: (a) Let i e S k . If i £ H k there is a point V E A k such that |v| = d (A k ) and v,- ^ 0. Then 
|v'| < d (Ajt). Assume now / ^ //^ and choose a point veAj such that |v| = d (Ak)- Then 
\v'\=d(A k ). 

(b) If i £ H k then there would be a point v G A^ with |v'| < d (A^) , which is in contradiction with 
the assumption d (A^) = d (A^) . Therefore i ^ Further for all points V EA k with |v| = d (A^) 
we have |v'| = J (A*) , therefore H k C 

(c) is a direct consequence of (48). 

(d) is verified by a direct calculation. 

Lemma 2 is the key to prove that the sequence (51) decreases: According to lemma 2(a), d (A k ) 
either decreases or remains constants. If it remains constant, lemma 2(b) states that (#4 — #H k ) 
either decreases or remains constant. If also this number remains constant, lemma 2(c) guar- 
antees T k = T k ' and lemma 2(d) allows us to descend in dimension and to consider the simpler 
problem for A^+i. 

It is an easy exercise to complete the proof and to show that if d (A r ) = the choice of T r forces 
d (A r ) to decrease. 
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